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Abstract 

A two-dimensional bi-layer, square lattice Heisenberg model with different 

intraplane(J||) and interplane( Jj_) couplings is investigated. The model is first 

solved in the Schwinger boson mean-field approximation. Then the solution 

is Gutzwiller projected to satisfy the local constraint that there should be 

only one boson at each site. For these wave functions, we perform variational 

Monte Carlo simulation up to 24x24x2 sites. It is shown that the Neel order is 

destroyed as the interplane coupling is increased. The obtained critical value, 

J±/J\\ = 3.51, is smaller than that by the mean-field theory. Excitation 

spectrum is calculated by a single mode approximation. It is shown that 

energy gap develops once the Neel order is destroyed. 
75.10.-b, 75.10.Jm, 75.50.Ee 
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I. INTRODUCTION 



Spin pseudogap observed in underdoped YBa 2 Cu 3 07_ :r is one of the fascinating char- 
acters of the high-T c cuprates. NMR experiments showed that even above the transition 
temperature of superconductivity, T c , static uniform susceptibility and the NMR relaxation 
rate,T 1; decrease with decreasing temperature.B Neutron scattering experiments showed the 
decrease of low energy magnetic excitation with decreasing temperature and found precursor 
of a finite spin gap.@ It has been pointed out that these astonishing experimental results can 
be explained provided that there is a spin pseudogap in the normal state of high T c mate- 
rials. These phenomena indicating the spin pseudogap, however, have not been observed in 
the L^-xSr^CuC^ systems.! Therefore it is speculated that the number of the CuC>2 layers 
between the insulating layers is essential for the formation of this gap, although a successful 
theory has not been presented.^ - ! 

It is conceivable that the finite concentration of holes affect the spin configuration and the 
excitation considerably. However, as a first step toward understanding of the spin pseudogap 
behavior, it is meaningful to study the properties of the bi-layer Cu02 system at zero doping. 
Namely we investigate a bi-layer square-lattice Heisenberg model of spin 1/2. 

H = J\\ a2 ^ i ^ a ' Si+w,a + J± 22 " ^*> 2 ' (1) 

i w o=l,2 i 

where w = x,y, i + w represents a site next to the site i in the w direction, and Si, a is a spin 
1/2 operator at site % in plane a. The nearest neighbor spins interact antiferromagnetically 
with intraplane coupling constant Jm and interplane coupling constant J±. What we want 
to know is how the properties of the system changes as J±/J\\ = a increases: at what value 
of a the Neel order is destroyed, and how the excitation spectrum varies. 

As for the zero-temperature critical value of a c for the destruction of the Neel order, 
there have been several investigations by various methods: spin wave approximation given 
by Matsuda and Hida§0 and the Schwinger boson mean-field theoryl resulted in quite large 
critical value, a c = 4.24 for the former and 4.48 for the latter. On the other hand, more 



2 



sophisticated methods have resulted in much lower critical value. Quantum Monte Carlo 
calculation gives it as 2.51 ± 0. 010 and the dimer expansions, which is an approach from 
the a — > oo limit, gives 2.56.0 One of the aim of the present paper is to obtain this critical 
value by another method, the Schwinger-boson Gutzwiller-projection method. 

In Schwinger-boson Gutzwiller-projection method we first solve the Hamiltonian by 
Schwinger-boson mean-field theory. The obtained ground-state wave function is Gutzwiller 
projected to fix the spin at each site to be 1/2. The wave function thus obtained is a kind 
of RVB00 wave function where long-range bonds are allowed with amplitude depending on 
the distance between the sites.HH This method was first used by Chen and Xiu for the square 
lattice antiferromagnetic Heisenberg model.0 It was shown that the wave function obtained 
this way is quite close to the true ground-state. This method has also been applied to the 
anisotropic Heisenberg model. There it was shown that even in the one- dimensional limit 
the ground-state energy, — 0.4377 J per site , is quite close to the exact value, — 0.4431J per 
site.0 Therefore we expect that this method gives wave functions quite close to the actual 
ground-state in the present system, too. A merit of the present method is that the wave 
function is given as an RVB wave function. Thus vertically coupled dimer state in the limit 
of a — > oo, disordered state in the intermediate value of a, and the Neel state at small a 
can be described in a unified way by wave functions with the same structure. 

In this paper, using this method we show that the Neel order at small a is destroyed at 
a c = 3.51. It is expected that gap appears in the excitation spectrum at a > a c . This is 
confirmed by our calculation of the spectrum by a single mode approximation. To obtain 
these results we solve the present Hamiltonian by the Schwinger-boson mean-field theory 
in Sec.[n[ The obtained ground-state wave function is Gutzwiller-projected in Sec.|HT[ The 
single mode approximation for the RVB wave function is discussed in Sec.|V[ In Sec.|V|, 
we perform variational Monte Carlo simulation for these wave functions and calculate its 
energy, spin-spin correlation, staggered magnetization and low-lying excitation spectrum. 
In Sec.|VT| critical point and the excitation spectrum are discussed. 
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II. MEAN-FIELD SOLUTION 

We introduce four kinds of bose operators, Sj i(Ii f , s iAi i{a = 1,2) to express the spin 
operators 

S£a = S i,a,1 S i,a,l > = 2 ( S i,a,T S *AT — s i,a,| S i,o,l)- (2) 

The commutation relations of the spin operators Si are satisfied in this replacement. We 
impose a constraint, 

in order to guarantee 5 = 1/2. Then Hamiltonian is rewritten as follows; 



i W 0=1,2 CT 



+ o ^-L ^ y^( s t.l.(7 S i.2.-(7 S t,2,q- s i,l,-g s i,l,<T S i,2,-<T S i,2,-o- s i,l,cr) 

+ S^AaW- ( 4 ) 

i 0=1,2 cr 

Here \i is a chemical potential introduced to enforce the constraint Eq.(^) on the average. 
To solve the Hamiltonian in the mean-field approximation, we introduce the following mean- 
field order parameters A w>a , A z , and n a ^ , which give the amplitudes of the intralayer singlet 
correlations, interlayer singlet correlations, and an averaged occupation number, respectively, 

A w = A Wj 2 = — Aw, l = 2 ( S «,2,| S «+«',2,T — Sip^Si+wp,]) , (5) 

A z = 2< s *,i.i s i,2,T ~ s i,i,T a *,2,4)' ( 6 ) 

n a,a = ( S i,a,cr S i,a,cr) = ~- (7) 

After decoupling the Hamiltonian, we rewrite the operator using its Fourier transformation: 

V A fe 

where N is the total number of lattice sites for each layer, and k summation is taken over 
the Brillouin zone —it < k x < 7r, — 7r < k„ < n. The mean-field Hamiltonian Hmv is written 
as follows 



H mf = E E A (4,a,T s M,T + sL fc>a>i s_ fc)0ii ) + «7fc(4,2,T s -fe,2,| ~ 4,i,T s -M,l) 

fc a 

- ^7fc(sfc,2,TS-fc,2,| - Sfc,i,T s -fe,i,l) - «^(sfc,i,t s -fc,2,| - Sfc,2,TS-fc,l,l) 

+ (4,i,T s -*Ai ~ s l,2,T s -fc,i,i) + const - ( 9 ) 



with 



A = /i-J||, (10) 

7 fe = 2J||(A x sin/i; x + A y sin k y ), (11) 

5 = -U±A Z . (12) 

The Hamiltonian can be diagonalized by a paraunitary Bogoliubov transformation 

s M)T = (cosh 0£afc T - cosh ^/3 feT + isinh^o;l H -isinh^/3i fci ), (13) 

s -fc,i,l = -^(^ inh ^fc «fcT - isinh^/3 feT - cosh0+a:L fei + cosh 0^1^), (14) 

Sk,2,i = —i= (cosh 9£akt + cosh O^Pki + isinh fl^at^ + i sinh 0^7 /3l fe jJ, (15) 
V2 

s -fc, 2 a = ^(-«sinh0+a feT -isinh^/3 feT + cosh 0+aL fei + cosh ^/3l fei ), (16) 

where 



sinh^ = -^^^sgn( 7fe±( 5), (17) 



£ fe± = ^-(Tfci^) 2 . (18) 
After the transformation, Hamiltonian finally becomes 

Hmf = Yl E k+(®U a W + al kl a- kl ) + £ fe _(/3j T /3 fcT + f3[_ H P- kl ) + const. (19) 
fc 

The ground-states \G) is defined as the vacuum of the Bose operator a k ^, oc- k i, (3 k ^, and 
P- ki , such that a kr \G) = a- kL \G) = (3 k ^\G) = p. kl \G) = 0. 
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For a finite-size system, the self consistent equations for X,A x ,A y , A z are given by Eqs.(j^- 
0), which lead to 

We find that the free energy takes the same minimal value for A x = A y (s-wave) and A x = 
— A y (d- wave) .@@ Since either state gives the same result, we consider only the s-wave state 
from now on. We denote A^ = A y = An and — iA z = Aj_. The solution depends on the 
size of the system N. When N is finite, E k ± never becomes zero. However, in the limit of 
iV — > oo it is possible that vanishes at k — K± = ±(tt/2,tt/2). In such a case it is 
known that we need to introduce the Bose condensate rig, and Eqs.( |20| - p2"D are rewritten as 

1 = w L L (w + + -t) dKdk ' + n » (23) 

sm k w — h — dk x dky + n B , (24) 



11 4(2tt) 2 7-.7-. ^ \ E k+ £ fe _ 

^ = 4^ L L " ^r) + nB - (25) 

When the Bose condensate ub becomes finite, we have A = 4J||A|i + J±A±. The self- 
consistent equations are numerically solved. Figure 1 shows the a dependence of order- 
parameters Ay, Aj_, Bose condensate ub, and energy gap E g . Bose condensate vanishes at 
a = 4.48, and the gap opens for a > 4.48. The intralayer RVB order parameter, An, vanishes 
at a = 4.62. For a > 4.62, only the interlayer nearest neighbor spin-spin correlation exists. 

The intra(inter) layer spin-spin correlation (Si, a • Sj, a )((Si,i ' Sj,2)) i n the grand-state 
are given as 
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1 ^, A A 



(Si, a • Sj, a ) = - [— + C ° S k ' n « + HB COS K + ' Vi 'i 



k 



2 



\ ^Tn ^^jT + ^k~^ sin k ' ri ' j + nB sin K+ ' ri J 2 ' ^ 



[Si,i • Sj, 2 ) - 



A , 2 

sin k ■ r itj + ra B sin K + ■ r^j 



3r 1 v-,7fc + <^ Ik-S, n2 
— > ( ) cos k • r. ,• + n r cos A. • r» g - 



(27) 



where r^j = — r^. The summations over k in Eqs. ( PB| , |2"7D vanish in the limit |rjj| — > oo. 
Therefore, the correlation extends to infinity only if ub > 0, which means the existence 
of antiferromagnetic long-range order. Thus, in the mean-field approximation, the critical 
point of order-disorder transition is 4.48. 

III. RVB WAVE FUNCTION 

The ground-state wave function obtained in the mean-field theory is expressed as 



\G) = Y[exp 

k 



-.tanh#^ + tanh# fc ^ j ^ ^ N 

1 o \ S k,2,-\ S -k,%[ ~ S k,l,\ S -k,l,lJ 



.tanh^ — tanh# fe . j j- ; t . , 



(28) 



where |0) is the vacuum of the Schwinger bosons. By the Fourier transformation for 
s\ip s -k,i,i> s i,2,T> s -fc,2,!' we can S e ^ a real-space representation for this ground-state, 



\G) = exp [J2 otj (4,2,T s i,2a - 4i,t4i,i) + M4u s Ut ~ 4i,t s ki)] 1°) 



(29) 



— 53[tanh^ + tanh0j 



2A 



Y^[tanh6£ - taring 



exp(ifc • rij), 
exp(ifc • r i;i ). 



(30) 
(31) 



It is evident that the local constraint, Eq.(||) is not satisfied in this wave function. We remove 
this difficulty by projecting the wave function to a space where each site is singly occupied. 
Namely, we perform the Gutzwiller projection, using Gutzwiller projection operator P, 



N . 



|0>. 



(32) 



From Eq . (p2|) it is clear that the ground-state \G) is an RVB state. The weights of the bond, 



ciij and bij, decay proportional to r i 3 except for b it j at small J±. Although it would be 



possible to regard every ay and &y as variational parameters, we here restrict them to be 
those given in Eqs.(^,^l|). In the case of a = 0, this restriction is justified by the result 
itself: Chen and Xiu@ have shown that this choice of ay gives excellent results for the 
ground-state energy and the staggered magnetization. The weights ay and fey depend on 
a = Jj_/ J\\ through the order parameters. We consider this a in ay and fry as a variational 
parameter. In order to avoid confusion, we use a symbol a p to mean the value of a used to 
obtain the ground-state. 



IV. EXCITATION SPECTRUM 

Once the approximate ground-state is obtained, excitation spectrum can be calculated 
by a method given by Feynman for liquid 4 He, namely the single mode approximation.illil 
The essential point of this method is to consider a low-lying excited state intuitively and 
calculate excitation spectrum from a known ground-state. In our case, low-lying state of 
this Hamiltonian should be the spin wave excitation. Thus we consider the following excited 
states. 

\E ± ) = (S^±S k ^)\G) , (33) 

^k,a = ~7Ji^ E ^i,a e% ^ T% i (34) 



where \E±) is variational excited states. Excitation spectrum, u±(k), is calculated as 

- ± (*) = t! • (35) 

S ±« = I E<G|(SJ ± S+ )(S-! ± S- 2 )\G)e ,k r '.i , (36) 
f±( k ) = ^ E(G\(Si ± S&[H, (Sr, ± S- i2 )]\G)e ik - r ij 



N 



iT^^E^K^ii^^i^-^^i^e^ . (37) 



N , 



Here, u = ±x, ±y, i + to' represents a site next to the site i in the <J direction, S±(k) is 
the static structure factor and f±(k) is a 3-point correlation function of spin operators. The 
two modes represent in-phase, u + (k), and out-of-phase, u-(k), spin excitations of the two 
layers. 

Since \G) is an RVB state, we must consider a loop covering associated with two valence 
bond configurations, |ci), |c 2 ), to calculate Eqs. (|36|j37|) .El For Eq.(|36|) we use known results, 

- (i, a), (J, b) belong to the same loop and the same sub-lattice. 



{ C l\Si,aS jyb \c 2 ) 

(ci\c 2 ) 



(i, a), (j, b) belong to the same loop and different sub-lattice. 



(i, a), (j, b) belong to the different loop. 



(38) 



For Eq . (pTj) the following rule is found, 

1 



\ c l I S% , b Sj+ 8, c I c 2 ) 

(ci|c 2 ) 



[i, b),(i + 5, c) belong to the same loop 



1 
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and (I, a) = (i + 5, c). 

(i, b),(i + 5, c) belong to the same loop (39) 
and (I, a) = (i, b). 
otherwise. 

Here, i + 5 means the nearest neighbor of i-th site. S±(k) can be calculated directly from 
the first rule. Using the second rule, f±(k) becomes 



f±(k) = ^(2-cosk x ~cosk y )J2i: E (G\S+ +w , a Sr, a S? +w JG) 

iV i w a=l,2 



(40) 



Thus we have only to count the number of the nearest neighbors in the same loop for each 
loop covering. This simplifies the numerical calculation. 



V. NUMERICAL RESULTS 



In this section, we show numerical results of the ground-state energy, spin-spin correla- 
tion, staggered magnetization, and the excitation spectrum as a function of a. We perform 
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Monte Carlo simulations in which RVB states are sampled to satisfy detailed balance for 
lattices with LxLx2 sites, where L < 24. All the numerical calculations are performed with 
periodic boundary conditions. For each system size we solve the self-consistent equations 
( p0| - |2~2| ) , and calculate a^, bij to be used to construct the wave function at that system size. 

A. Ground-state energy 

The energy per site of the bi-layer Heisenberg model, E, is given by the nearest-neighbor 
spin-spin correlations €\\(L,a p ) and e±(L,a p ) for a given wave function specified by the 
parameter a p . 

E(L,a p ) = 2J||e||(L,a p ) + ^J±e ± (L,a p ) , (41) 

where 

^^) = 775EEE( G I S m.-WG) , (42) 

4iv i w a=l,2 

e^a p ) = ±- 2 Y.(G\S hl -S h2 \G) . (43) 

To estimate the energy in the thermodynamic limit, the size dependence is examined 
and we find the following size scaling for any fixed a p , 

e ll (L,a p ) = e ll (a p )+XL- 3 + --- , (44) 

e±(L,a p ) = e ± (a p ) + XL~ 3 + ■ ■ ■ , (45) 

where A is a constant. This size-scaling coincides with the spin wave theory for a square 
lattice. In Fig. 2, e\\(a p ) and e±(a p ) are shown. Open circles and solid circles indicate 
en (o!p) and ej_(a p ). Error bars show the standard deviation of the Monte Carlo simulation. 
The interplane nearest-neighbor spin-spin correlation en has a value of —0.3333 ± 0.0006 at 
a p = 0, which is quite close to the best estimated value of — 0. 3348.1111 The magnitude 
of ey decreases as a p increases and finally vanishes at a p = 4.62. On the other hand, 

10 



the magnitude of e±(a p ) increases and saturates to 0.75 at a p = 4.62. At a p > 4.62 the 
intraplane spin correlation vanishes and dimerized state is realized. 

The ground-state energy per site at a given a is calculated as a minimum of E(a p ) = 
2J||e||(o!p) + | J±e±(a p ) with respect to a p . Thus, we can get an optimal variational parameter 
and energy for a given a. The relation between the variational parameter(o;p) and a real 
coupling(o;) is shown in Fig. 3, and the ground-state energy per site is shown in Fig. 4. In 
Fig. 4 we also show the energy of dimerized state per site, — |a Jy, (straight line) for reference. 
The difference between the optimal energy and the dimerized energy becomes smaller with 
increasing interlayer coupling. 

B. Staggered magnetization 

We calculated the spin-spin correlation, (Si >a • Sj^), between arbitrary two sites (i,a) 
and (j,b). The results for 24 x 24 x 2 lattice system are shown in Fig. 5, where absolute 
value of the intralayer spin-spin correlation is plotted as a function of the distance between 
the two sites. Open circles are for a = 0.4 and solid circles are for a = 4.6. It is obvious that 
there is a long-range order at a = 0.4 and no long-range order at a — 4.6. In the latter case, 
the correlation decreases exponentially and the typical correlation length for the disordered 
state is of the order of a lattice constant. 

The long-range order is of the antiferromagnetic one. In the ordered phase staggered 
magnetization of the infinite size system is obtained from size dependence of the staggered 
spin-spin correlation between the mostly separated pairs. For a given lattice size L, we 
calculated both the intralayer correlation M (L) 2 , and interlayer correlation Mi(L) 2 : 

M o(Lf^^-Z E(\Si,a-Sj, a \) , (46) 

ZiV i,j a=l,2 

mi(l)* = ±y;(\ s v s *\) > ( 47 ) 

where the summation is taken for all the pair of % and j such that Ti — rj = (±L/2, ±L/2). 
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Except at a = 0, Mq(L) and M±(L) coincide within the Monte Carlo statistical error. 
As shown in Fig. 6, they are well fitted by the size scaling, 

M (L) = M X {L) = M(oo) + ■ •• , (48) 

where fi is a constant. This scaling agrees with the prediction of the spin wave theory and 
arguments given by Huse.il The staggered magnetization Mq = M(oo) as a function of a is 
given in Fig. 7. In this figure, the results of the mean-field theory(MFT) are also shown. In 
the case of small a, the interlayer coupling enhances the antiferromagnetic long-range order. 
This is because the system acquires a weak three-dimensionality and quantum fluctuation 
is suppressed. On the other hand, for larger a, the magnetizations are suppressed. This 
behavior is consistent with the result of Matsuda and Hida in the spin-wave theory! The 
staggered magnetization vanishes at a c = 3.51 ± 0.05. 

C. Excitation spectrum 

We calculate the structure factor, S±(k), and excitation spectrum, u±(k) as a function 
of coupling a. The calculations are done for 24 x 24 x 2 lattice. The behavior of S±(k) 
and u}±(k) strongly depends on whether the system has long-range order or not. The result 
of S±(k) is shown in Fig. 8 and u±(k) is shown in Fig. 9. Here, three typical couplings 
are taken; a = 0.4 (open circles), a = 2. 4 (closed circles), a = 3. 6 (open squares). The third 
coupling is for the system in the disordered phase. For each figure, (a) is for the plus mode 
and (b) is for the minus mode, and T = (0, 0), X = (0, it), M = (ir, 7r) in momentum space. 

It is obvious from Fig. 8 that S + (k) of the ordered state(a = 0.4, 2.4) is proportional to 
k near T point and S-(k) has an antiferromagnetic peak at M point. On the other hand, 
S+(k) at a = 3.6 increases quadratically with k near T point. (See inset of Fig. 8(a)). In the 
Neel state, the excitation is gapless at two points. One is w + (k) at F point. Around this 
point, since the function f+(k) in Eq.([$5|) behaves like f+(k) oc k 2 and the structure factor 
does S + (k) oc k, the excitation is proportional to k. The other is a>_(fe) at M point where 
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S-(k) diverges due to the antiferromagnetic long-range order. Thus, the gap opens when 
the structure factor becomes proportional to the square of k for the former point and when 
the structure factor does not diverge, that is, the system becomes the disordered state for 
the latter point. In the former case, we should determine the critical coupling, a c2 , where 
the gap opens. We take five k x points and do the following fitting along the T — X line; 
S+(k x ) = a\k x + + a^k^, where a±, 02, and a 3 are fitting parameters. The result of a± 
versus a is shown in Fig. 7 (open squares). Comparing the result of staggered magnetization 
with this coefficient, we find that the critical point a C 2 is equal to the a c within the statistical 
and fitting errors, which gives a c = 3.51 ± 0.05. The a dependence of the gap is shown in 
Fig. 10. All values are scaled by Jy. Open circles are for c<j + (0,0) and closed circles are for 
c<j_(7r,7r). In the disordered phase, excitation energy c<j_(7r,7r) always takes smaller value. 
Spin wave velocity along the T — X line is calculated for the ordered state and the result 
is shown in the inset of Fig. 10. Here, Z c is the renormalization factor. Namely spin wave 
velocity is given by \plZ c J\\. As the coupling increases, the velocity first slightly decreases 
and then suddenly increases near the critical point. 



VI. DISCUSSIONS 

In this paper we first solved the Hamiltonian by the Schwinger-boson mean-field theory. 
Then the solution is Gutzwiller projected to obtain variational ground state wave functions, 
which are examined by the Monte Carlo simulation for finite sizes. 

We first see the advantage of our variational Monte Carlo simulation. In the mean-field 
calculation, the system becomes dimerized for a > 4.62. For this region, interlayer order 
parameter, Ay, is zero and only dimer coupling between the layers is permitted. In addition, 
excitation spectrum becomes flat in momentum space; E g {k) = y/X 2 — 5 2 . As a matter of 
fact, theoretically, this state must be only realized at a — > 00. This disadvantage is removed 
in Monte Carlo simulation. It is estimated from Fig. 4 that the virtual critical point where the 
system stabilizes with the dimerized state is 11.0. This means that the Gutzwiller projection 
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improves the wave functions. The improved wave function can describe the disordered state 
without dimerization at least up to a = 11.0. 

There have been many investigations for the order-disorder critical point. Our mean-field 
result is essentially same as the previous report^ and the modified spin wave theory.0 These 
give the critical value of a around 4.5.0 This value is much larger than the results by the 
other methods: 2.56 by the dimer expansion, and 2.51 ± 0.01 by the quantum Monte Carlo 
method. However, these latter values are still formidably larger than the value of a realized 
in the bi-layer cuprates. Our motivation for this work was to see if our method gives the 
critical value closer to the experimental value or not. Our result, a c = 3.51 ± 0.05,0 is 
not for this expectation, and confirms the previous theories that without doping bi-layer 
Heisenberg model will not give an explanation for the spin gap behavior of the experiments. 

We also calculated the excitation spectrum, especially for the disordered phase. It is not 
obvious whether the system has always a finite gap in disordered state. For instance, there 
is no long-range order for the one-dimensional s — 1/2 antiferromagnetic Heisenberg model, 
though the excitation spectrum is gapless. In our bi-layer two-dimensional Heisenberg model, 
we find that there is always a finite gap for the disordered state. Within the statistical and 
fitting errors, it occurs at a c = 3.51 ± 0.05. In disordered region, the spin-spin correlation 
decays as the distance exponentially. The structure factor, S-(k), near the critical coupling, 
however, has a large maximal value at M point, which makes the excitation spectrum be 
minimized at that point. This shows that even in the disordered state the antiferromagnetic 
spin fluctuation is strong. It should be remarked that even though the spectrum u_(k) at 
a = 3.6 looks singular at Q = (tt, it), this is not the case. Around Q it should be quadratic 
in {k — Q). Such a behavior is not apparent in Fig. 9(b) due to the lack of data close enough 
to Q. 

At a = where the model becomes the single layer Heisenberg model, our result can 
be compared with other methods: spin wave theory, series expansions and single mode 
approximation.&B Our result of the spectrum is roughly proportional to those of other 
methods over the entire Brillouin zone. The maximal value is around 2.65Jy at X or L = 
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(7r/2, 7r/2) point. Series expansions predicted the maximum is about 2.35J|| at L pointed 
and single mode approximation based on the expansions around the Ising limit estimated 
the maximum about 2.5Jii at L point.0 Both results are close to our result. The most 
remarkable difference from the other methods is the spin wave velocity. The renormalization 
factor, Z c , is 1.99 ± 0.03 at a = which is 1.69 times larger than the best estimated value 
around 1.18±O.O2.0 This difference indicates that multi-magnon contribution to S+(k) is not 
negligible. However, since this method gives qualitatively correct behavior, we believe it gives 
qualitatively correct spectrum at a > also. Finally we remark that the non-monotonous 
behavior of spin wave velocity with increasing interlayer coupling can be understood from 
that of the coefficient (g^ ) of the structure factor shown in Fig. 7, since the spin wave velocity 
is inversely proportional to a\. 

In conclusion, we have investigated the bi-layer Heisenberg model using the Schwinger- 
boson Gutzwiller-projection method. We find that there is an order- disorder transition with 
increasing interlayer coupling. The critical point is a c = 3.51 ± 0.05. Excitation spectrum 
can be calculated for wide range of coupling and we find that the spin excitation has always 
a finite gap for disordered phase and the minimum of the spectrum is located at M point. 
Our model corresponds to the half-filled case for high-T c cuprates. Although a c in this case 
is quite large, it is possible that hole doping reduces the value extremely. Then it will be 
possible that our disordered state continuously changes into the spin gap state. The similar 
treatment for a hole doped model, t — t' — J model, is our next problem. 
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FIGURE CAPTIONS 



FIG. 1. Mean-field values of order parameters Ay, Aj_, Bose-condensate n B and energy 
gap E g as a function of a. 

FIG. 2. The nearest neighbor spin correlation for each direction is shown. Open circles 
are for ey, and solid circles are for e±. Error bars result from Monte Carlo statistical errors. 

FIG. 3. The variational parameter a p which minimizes the ground state energy for a 
given parameter a. 

FIG. 4. Total energy per site as a function of a. Open circles are for variational Monte 
Carlo results and straight line is for dimerized state, —0.375a. 

FIG. 5. Spin-spin correlation for a = 0.4(open circles) and a = 4.6(solid circles). Each 
calculation is done for 24 x 24 x 2 lattice. Here, r it j means the distance between two sites. It 
is obvious that there is a long-range order for a = 0.4 but no long-range order for a = 3.6. 

FIG. 6. M (L) versus l/L for a = 0.0, 0.8, 1.7, 3.1, and 4.6. 

FIG. 7. Staggered magnetization as a function of a. Open circles are for the mean- 
field theory, and solid circles are for variational Monte Carlo results. The magnitude of 
the /c-linear term in the expansion of S+(k) around the T point, ai, is also shown by open 
squares. 

FIG. 8. The structure factors (a)S + (k) and (b)S-(k). Open circles, closed circles, and 
open squares are for a = 0.4, 2.4, and 3.6, respectively. Inset shows the detailed structure 
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of a = 3.6 along the T — X line. Note that the value at M point of S-(k) is too large to be 
shown in the figure. 

FIG. 9. Excitation spectrum (a)c<j + (fc) and (b)a>_(fe). The same values for a are chosen 
and indicated by the same symbols as in Fig. 8. At a — 3.6, gap opens at T point for u>+(k) 
and at M point for ou_(k). 

FIG. 10. The a dependence of the gap for cj + (0,0) and c<j_(7r,7r). In the inset, the 
renormalization factor of the linear spin wave velocity, Z c , is also shown. 
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